rm(list=ls())
# setwd(...) # set to appropriate working directory
data <- read.csv(file="MainData.csv")

# install necessary packages: install.packages("package name")
library(aod)
library(sandwich)
library(stargazer)
library(lmtest)
library(MASS)
library(ordinal)
library(multiwayvcov)
library(car)
library(rms)
library(lme4)

#########################################################################################################################
### Table 1 ###
# Model 1: log(Casualty)
dummy.matrix <- as.matrix(cbind(as.numeric(data$Year == 2006), as.numeric(data$Year==2010), as.numeric(data$Year==2014)))

Log_Casualty.Model <- lm(Ethnic_Vote_Share ~ Log_Casualty:dummy.matrix+dummy.matrix+
                   Municipality-1, data)
vcovCL1 <- cluster.vcov(Log_Casualty.Model, ~data$Municipality+data$Year)
coeftest(Log_Casualty.Model, vcovCL1) # results
nobs(Log_Casualty.Model) # 428 observations
summary(Log_Casualty.Model)$r.squared # R-squared: 0.9905888

# Model 2: Refugees
Refugees.Model <- lm(Ethnic_Vote_Share ~ Refugees:dummy.matrix + Municipality
                     + dummy.matrix-1, data)
vcovCL2 <- cluster.vcov(Refugees.Model, ~data$Municipality+data$Year)
coeftest(Refugees.Model, vcovCL2) # results
nobs(Refugees.Model) # 428 observations
summary(Refugees.Model)$r.squared # R-squared: 0.9908892

# Model 3: Prison Camps
Prison_Camps.Model <- lm(Ethnic_Vote_Share ~ Prison_Camps:dummy.matrix
                         + dummy.matrix+Municipality-1, data)
vcovCL3 <- cluster.vcov(Prison_Camps.Model, ~data$Municipality+data$Year)
coeftest(Prison_Camps.Model, vcovCL3) # results
nobs(Prison_Camps.Model) # 388 observations
summary(Prison_Camps.Model)$r.squared # R-squared: 0.9906923

# Model 4: Claims 
Claims.Model <- lm(Ethnic_Vote_Share ~ Claims:dummy.matrix
                   + dummy.matrix+Municipality-1, data)
vcovCL4 <- cluster.vcov(Claims.Model, ~data$Municipality+data$Year)
coeftest(Claims.Model, vcovCL4) # results
nobs(Claims.Model) # 396 observations
summary(Claims.Model)$r.squared # R-squared: 0.9904656
#########################################################################################################################

#########################################################################################################################
### Appearance of Social Democratic Party of BiH (SDP BiH) on ballot in 2016 local elections (footnote 31) ###
SDP_Appearance <- read.csv("SDP_Appearance.csv")
SDP_Appearance$SDP_Not_Appear_Pop <- ifelse(SDP_Appearance$SDP_Appear == 1, 0, SDP_Appearance$Population)
(sum(SDP_Appearance$SDP_Not_Appear_Pop)/sum(SDP_Appearance$Population))*100 # 12.29336
# The SDP BiH did not appear on the ballot in municipalities that constitute roughly 12% of the population.
#########################################################################################################################

#########################################################################################################################
### F-Tests for difference between 2006 and 2010, 2010 and 2014, 2006 and 2014 (footnote 41) ###
linearHypothesis(Log_Casualty.Model, c('Log_Casualty:dummy.matrix1=Log_Casualty:dummy.matrix2'), vcov.=vcovCL1) # 2006 and 2010
linearHypothesis(Log_Casualty.Model, c('Log_Casualty:dummy.matrix2=Log_Casualty:dummy.matrix3'), vcov.=vcovCL1) # 2010 and 2014
linearHypothesis(Log_Casualty.Model, c('Log_Casualty:dummy.matrix1=Log_Casualty:dummy.matrix3'), vcov.=vcovCL1) # 2006 and 2014
# Results: 2006 and 2010 not distinguishable, 2010 and 2014 distinguishable, 2006 and 2014 distinguishable
#########################################################################################################################

#########################################################################################################################
### Substantive Effects of Violence Measures (in text) ###
# log(Casaulty) - Table 1, Model 1
diff(range(data$Log_Casualty)*3.13527) # difference from lowest to highest in 2006
diff(range(data$Log_Casualty)*2.06516) # difference from lowest to highest in 2010
diff(range(data$Log_Casualty)*0.41477) # difference from lowest to highest in 2014
quantile(data$Casualty, probs = c(.25, .75)) # 25% quantile: 1.125952, 75% quantile: 2.345260
diff(quantile(data$Log_Casualty, probs = c(.25, .75))*3.13527) # difference from 25% to 75% in 2006
diff(quantile(data$Log_Casualty, probs = c(.25, .75))*2.06516) # difference from 25% to 75% in 2010
diff(quantile(data$Log_Casualty, probs = c(.25, .75))*0.41477) # difference from 25% to 75% in 2014

# Refugees - Table 1, Model 2
quantile(data$Refugees, probs = c(.25, .75)) # 25% quantile: 1.620941, 75% quantile: 5.157595
diff(quantile(data$Refugees, probs = c(.25, .75))*1.64919) # difference from 25% to 75% in 2006
diff(quantile(data$Refugees, probs = c(.25, .75))*1.49478) # difference from 25% to 75% in 2010
diff(quantile(data$Refugees, probs = c(.25, .75))*1.18340) # difference from 25% to 75% in 2014

# Refugees - Table 1, Model 3
quantile(data$Prison_Camps, probs = c(.25, .75), na.rm=TRUE) # 25% quantile: 0.06222775, 75% quantile: 0.24421350
diff(quantile(data$Prison_Camps, probs = c(.25, .75),na.rm=TRUE)*13.03477) # difference from 25% to 75% in 2006
diff(quantile(data$Prison_Camps, probs = c(.25, .75),na.rm=TRUE)*13.51427) # difference from 25% to 75% in 2010
diff(quantile(data$Prison_Camps, probs = c(.25, .75),na.rm=TRUE)*10.86287) # difference from 25% to 75% in 2014

# Claims - Table 1, Model 4
quantile(data$Claims, probs = c(.25, .75), na.rm=TRUE) # 25% quantile: 2.056973, 75% quantile: 5.770101
diff(quantile(data$Claims, probs = c(.25, .75),na.rm=TRUE)*1.27168) # difference from 25% to 75% in 2006
diff(quantile(data$Claims, probs = c(.25, .75),na.rm=TRUE)*1.52613) # difference from 25% to 75% in 2010
diff(quantile(data$Claims, probs = c(.25, .75),na.rm=TRUE)*1.06479) # difference from 25% to 75% in 2014
#########################################################################################################################

#########################################################################################################################
### In-Municipality and Out-of-Municipality Voting ###
# 1997 Election
election1997 <- read.csv("Election1997.csv")
(sum(election1997$HDZ_In, election1997$SDA_In, election1997$SDS_In)/sum(election1997$Total_In_Votes))*100
# 59% of in-municipality votes received by SDA, HDZ BiH, or SDS in 1997.
(sum(election1997$HDZ_Out, election1997$SDA_Out, election1997$SDS_Out)/sum(election1997$Total_Out_Votes))*100
# 83% of in-municipality votes received by SDA, HDZ BiH, or SDS in 1997.

# 1998 Election
election1998 <- read.csv("Election1998.csv")
election1998$In_Votes <- election1998$In_Votes_FBiH + election1998$In_Votes_RS 
election1998$Out_Votes <- election1998$Out_Votes_FBiH + election1998$Out_Votes_RS 
(sum(election1998[election1998$Party == "SDA" | 
             election1998$Party == "HDZ BiH" | 
             election1998$Party == "SDS", "In_Votes"])/sum(election1998$In_Votes))*100
# 51% of in-municipality votes received by SDA, HDZ BiH, or SDS in 1998.
(sum(election1998[election1998$Party == "SDA" | 
                    election1998$Party == "HDZ BiH" | 
                    election1998$Party == "SDS", "Out_Votes"])/sum(election1998$Out_Votes))*100
# 69% of out-of-municipality votes received by SDA, HDZ BiH, or SDS in 1998.
#########################################################################################################################

#########################################################################################################################
### Correlations between ethnic homogeneity and violence (footnote 59) ###
cor(data$HHI.1991, data$Casualty) # -0.1114354
cor(data$HHI_SB, data$Casualty)   # 0.14654
cor(data$HHI_CB, data$Casualty)   # 0.1177229
cor(data$HHI_SC, data$Casualty)   # -0.2097433
#########################################################################################################################

#########################################################################################################################
### Survey data analysis ###
SurveyData <- read.csv("SurveyData.csv")

Ethnic_Friends_Data <- na.omit(SurveyData[,c('Ethnic_Friends','Log_Casualty','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                             'Pop_Density','Income_PC','Municipality')])

Closest_Friends_Data <- na.omit(SurveyData[,c('Closest_Friends','Log_Casualty','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                              'Pop_Density','Income_PC','Municipality')])

National_Trust_Data <- na.omit(SurveyData[,c('National_Trust','Log_Casualty','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                              'Pop_Density','Income_PC','Municipality')])

Representation_Data <- na.omit(SurveyData[,c('Representation','Log_Casualty','Gender','Age','HHI_SB','HHI_CB','HHI_SC',
                                             'Pop_Density','Income_PC','Municipality')])

### Table 2 ###
# Model 1: Ethnic Friends
Ethnic_Friends.model <- lrm(Ethnic_Friends ~ Log_Casualty
                             + Gender
                             + Age
                             + HHI_SB
                             + HHI_CB
                             + HHI_SC
                             + Pop_Density
                             + Income_PC, y=TRUE, x=TRUE, data=Ethnic_Friends_Data)
robcov(Ethnic_Friends.model,cluster=Ethnic_Friends_Data$Municipality)

# Model 2: Closest Friends
Closest_Friends.model <- lrm(Closest_Friends ~ Log_Casualty
                            + Gender
                            + Age
                            + HHI_SB
                            + HHI_CB
                            + HHI_SC
                            + Pop_Density
                            + Income_PC, y=TRUE, x=TRUE, data=Closest_Friends_Data)
robcov(Closest_Friends.model,cluster=Closest_Friends_Data$Municipality)

# Model 3: National Trust 
National_Trust.model <- lrm(National_Trust ~ Log_Casualty
                             + Gender
                             + Age
                             + HHI_SB
                             + HHI_CB
                             + HHI_SC
                             + Pop_Density
                             + Income_PC, y=TRUE, x=TRUE, data=National_Trust_Data)
robcov(National_Trust.model,cluster=National_Trust_Data$Municipality)

# Model 4: Representation 
Representation.model <- lrm(Representation ~ Log_Casualty
                            + Gender
                            + Age
                            + HHI_SB
                            + HHI_CB
                            + HHI_SC
                            + Pop_Density
                            + Income_PC, y=TRUE, x=TRUE, data=Representation_Data)
robcov(Representation.model,cluster=Representation_Data$Municipality)
#########################################################################################################################


